Time Series Covid Project
Reagan Meagher and Megan Riley
#Modeling Covid19 Time Series
Covid19 is a worldwide pandemic that will likely define 2020. In the United States, currently over four million people have been infected and over 150,000 have died as a result of Covid19. As the pandemic continues, limiting infections, serious harm, and death is a primary concern for all involved.
As this is a novel illness, we know relatively little, but understanding how Covid19 is spreading and judging the severity of an outbreak can be approximated with the data we have available. In this report we aim to build effective time series models to forecast future Covid19 cases using the techniques we have learned from this Time Series course. ##Goal One: Data Collection
The data source we are using is sourced from The Covid Tracking Project.
initial_data_fl <- read.csv(file="https://raw.githubusercontent.com/megnn/TimeSeries_Covid/master/covidtracking_FL_data.csv", header=TRUE)
initial_data_us <- read.csv(file="https://raw.githubusercontent.com/megnn/TimeSeries_Covid/master/covidtracking_US_data.csv", header=TRUE)
initial_data_fl = initial_data_fl[order(nrow(initial_data_fl):1),]
initial_data_us = initial_data_us[order(nrow(initial_data_us):1),]
len_fl = dim(initial_data_fl)[1]
len_us = dim(initial_data_us)[1]Insert exploration of cases, deaths, hospitalizations etc.
Plotting the daily positive cases in Florida and US.
plot(x = seq(1,len_fl), y = initial_data_fl$positiveIncrease, type = "l")
plot(x = seq(1,len_us), y = initial_data_us$positiveIncrease, type = "l")
##Positive Percentage
Positive Percentage is a statistic that calculates daily positive tests as a percentage of daily overall tests returned. We calculated this column and added it to our data below followed by some visual exploration of the statistic itself.
Overall we see a clear instance of high and often 100% positive test rates early on in the first days and weeks of the pandemic spread. We understand this as a result of the fact that Covid19 spread fast and we had more community spread than anticipated early on without the testing available. It is abundantly clear that when we have extremely high percent positive rates near 100% we can expect true positive case numbers at the time to be under represented. But without better epidemlogical understanding we can’t make judgement calls on true case numbers when percent positives rise from 5% to 10% as we see begin to happen somewhat in recent days in Florida.
for (i in 1:nrow(initial_data_fl)) {
n <- round((initial_data_fl$positiveIncrease / initial_data_fl$totalTestResultsIncrease) * 100, digits = 4)
initial_data_fl$positive_percentage <- n
}
for (i in 1:nrow(initial_data_us)) {
n <- round((initial_data_us$positiveIncrease / initial_data_us$totalTestResultsIncrease) * 100, digits = 4)
initial_data_us$positive_percentage <- n
}
#Percent Positive Exploration
plot(x = seq(1:len_fl), y = initial_data_fl$positive_percentage, type = "l")
plot(x = seq(21:135), y = initial_data_fl[21:135,42], type = "l")
plot(x = seq(1:len_us), y = initial_data_us$positive_percentage, type = "l")
plot(x = seq(45:177), y = initial_data_us[45:177,26], type = "l")
Positive Percentage as a metric is a measure of two main things, how many tests are we administering and how many positives are we receiving. If tests are skyrocketing while positive cases are increasing, we would see a stable or even diminishing line which could indicate not a pandemic under control but simply better testing resources but could be interpreted as a pandemic managed.
Keeping tests increasing to continue to keep percent positives level is a good indication we have leveled up our resources to continue to diagnose the pandemic at the same level, but if we need to scale up our testing to keep the same positive percentage, there is more covid spread.
However, an increasing positive percentage is a good indicator that our testing resources may not be up to actually up to tracking the current stage of the pandemic.
plot(x = seq(1:len_us), y = initial_data_us$totalTestResultsIncrease, type = "l")
plot(x = seq(1:len_fl), y = initial_data_fl$totalTestResultsIncrease, type = "l")
####Data Preperation
In order to model new case numbers by day we set up dataframes with only our date and positive increase amount per day.
newcases_fl <- dplyr::select(initial_data_fl, c("date", "positiveIncrease"))
newcases_us <- dplyr::select(initial_data_us, c("date", "positiveIncrease"))####Checking for NAs
We can see with the missing value analysis below that we have no NAs present in our new case data.
#Checking for NAs
md.pattern(newcases_fl)## /\ /\
## { `---' }
## { O O }
## ==> V <== No need for mice. This data set is completely observed.
## \ \|/ /
## `-----'
## date positiveIncrease
## 135 1 1 0
## 0 0 0# No NAs present
#Checking for NAs
md.pattern(newcases_us)## /\ /\
## { `---' }
## { O O }
## ==> V <== No need for mice. This data set is completely observed.
## \ \|/ /
## `-----'
## date positiveIncrease
## 177 1 1 0
## 0 0 0#No NAs present###Florida Daily Cases:
plotts.sample.wge(newcases_fl$positiveIncrease)
## $autplt
## [1] 1.00000000 0.90351878 0.85112908 0.81870330 0.78426726 0.72368910
## [7] 0.70892248 0.68040351 0.66167307 0.59052464 0.55736019 0.54019962
## [13] 0.50346164 0.45903473 0.44018864 0.39175866 0.35750776 0.31032227
## [19] 0.29436806 0.25707296 0.21075421 0.16689463 0.15038715 0.11609674
## [25] 0.09827282 0.08510643
##
## $freq
## [1] 0.007407407 0.014814815 0.022222222 0.029629630 0.037037037
## [6] 0.044444444 0.051851852 0.059259259 0.066666667 0.074074074
## [11] 0.081481481 0.088888889 0.096296296 0.103703704 0.111111111
## [16] 0.118518519 0.125925926 0.133333333 0.140740741 0.148148148
## [21] 0.155555556 0.162962963 0.170370370 0.177777778 0.185185185
## [26] 0.192592593 0.200000000 0.207407407 0.214814815 0.222222222
## [31] 0.229629630 0.237037037 0.244444444 0.251851852 0.259259259
## [36] 0.266666667 0.274074074 0.281481481 0.288888889 0.296296296
## [41] 0.303703704 0.311111111 0.318518519 0.325925926 0.333333333
## [46] 0.340740741 0.348148148 0.355555556 0.362962963 0.370370370
## [51] 0.377777778 0.385185185 0.392592593 0.400000000 0.407407407
## [56] 0.414814815 0.422222222 0.429629630 0.437037037 0.444444444
## [61] 0.451851852 0.459259259 0.466666667 0.474074074 0.481481481
## [66] 0.488888889 0.496296296
##
## $db
## [1] 13.9443470 12.1738932 9.4960087 5.0020029 1.2010108
## [6] -1.0130101 -0.1506086 0.2039391 -0.3463438 -0.9726125
## [11] -3.3208756 -5.1446096 -4.6985377 -4.0076175 -3.0331597
## [16] -4.4568722 -6.5691185 -9.5735201 -7.8085456 -4.7203532
## [21] -5.1362477 -6.5033145 -6.5841302 -11.7008911 -15.7768594
## [26] -20.3454088 -14.7914168 -13.3311995 -8.6134078 -11.8506844
## [31] -9.4294479 -10.6796170 -10.8717772 -6.8765678 -5.9591468
## [36] -4.1063148 -6.4365181 -5.2994271 -8.7718362 -14.6806581
## [41] -10.7718204 -10.0417172 -11.3897236 -10.1069047 -12.8544662
## [46] -10.2709895 -10.2140102 -8.6959976 -7.5565160 -7.8110494
## [51] -10.8377249 -12.1847651 -18.4078010 -14.4929054 -22.1053741
## [56] -21.8671829 -18.6457477 -17.1312143 -15.2219426 -8.8827236
## [61] -18.0198672 -11.1478112 -18.8429668 -14.8396876 -10.8187051
## [66] -6.0417196 -6.9676002
##
## $dbz
## [1] 10.8952916 10.4208269 9.6290499 8.5210982 7.1053841
## [6] 5.4084900 3.4951614 1.4961088 -0.3815372 -1.9206026
## [11] -3.0320075 -3.7937247 -4.3311799 -4.7225622 -5.0088115
## [16] -5.2333565 -5.4517814 -5.7160206 -6.0608486 -6.5041876
## [21] -7.0550944 -7.7189474 -8.4927699 -9.3490346 -10.2124871
## [26] -10.9466178 -11.3848655 -11.4222535 -11.0879684 -10.5084533
## [31] -9.8212142 -9.1329494 -8.5201783 -8.0386115 -7.7282189
## [36] -7.6148214 -7.7102817 -8.0113599 -8.4961143 -9.1170132
## [41] -9.7927830 -10.4087161 -10.8436066 -11.0271353 -10.9840189
## [46] -10.8151075 -10.6393942 -10.5538767 -10.6249324 -10.8939499
## [51] -11.3823560 -12.0894051 -12.9791522 -13.9554605 -14.8397014
## [56] -15.4053463 -15.5118058 -15.2066656 -14.6399148 -13.9276436
## [61] -13.1243083 -12.2660250 -11.4056316 -10.6130044 -9.9570803
## [66] -9.4915816 -9.2507536###US Daily Cases:
plotts.sample.wge(newcases_us$positiveIncrease)
## $autplt
## [1] 1.0000000 0.9588275 0.9195777 0.8839579 0.8613438 0.8422725 0.8250617
## [8] 0.7941208 0.7575713 0.7092840 0.6736316 0.6484889 0.6360652 0.6158459
## [15] 0.5843954 0.5425539 0.4973999 0.4629802 0.4433251 0.4227875 0.4027823
## [22] 0.3756549 0.3451029 0.3099138 0.2830722 0.2684966
##
## $freq
## [1] 0.005649718 0.011299435 0.016949153 0.022598870 0.028248588
## [6] 0.033898305 0.039548023 0.045197740 0.050847458 0.056497175
## [11] 0.062146893 0.067796610 0.073446328 0.079096045 0.084745763
## [16] 0.090395480 0.096045198 0.101694915 0.107344633 0.112994350
## [21] 0.118644068 0.124293785 0.129943503 0.135593220 0.141242938
## [26] 0.146892655 0.152542373 0.158192090 0.163841808 0.169491525
## [31] 0.175141243 0.180790960 0.186440678 0.192090395 0.197740113
## [36] 0.203389831 0.209039548 0.214689266 0.220338983 0.225988701
## [41] 0.231638418 0.237288136 0.242937853 0.248587571 0.254237288
## [46] 0.259887006 0.265536723 0.271186441 0.276836158 0.282485876
## [51] 0.288135593 0.293785311 0.299435028 0.305084746 0.310734463
## [56] 0.316384181 0.322033898 0.327683616 0.333333333 0.338983051
## [61] 0.344632768 0.350282486 0.355932203 0.361581921 0.367231638
## [66] 0.372881356 0.378531073 0.384180791 0.389830508 0.395480226
## [71] 0.401129944 0.406779661 0.412429379 0.418079096 0.423728814
## [76] 0.429378531 0.435028249 0.440677966 0.446327684 0.451977401
## [81] 0.457627119 0.463276836 0.468926554 0.474576271 0.480225989
## [86] 0.485875706 0.491525424 0.497175141
##
## $db
## [1] 14.2491777 14.9471733 6.9591395 9.5150572 4.6664320
## [6] 1.2452886 2.2001263 -1.8153526 0.2487892 -1.7451730
## [11] -2.0281764 -2.9144478 -2.1941621 -5.2629156 -5.2055932
## [16] -5.7078490 -5.9605659 -8.5976595 -9.4776476 -9.1217112
## [21] -15.2098111 -11.0130098 -12.2778403 -8.7998267 -3.7259575
## [26] -0.2641817 -7.1652992 -4.1204548 -7.1212961 -11.2647439
## [31] -8.3267561 -13.1787605 -11.1231072 -14.6051895 -11.4235887
## [36] -13.9657569 -10.8725668 -12.1561887 -10.1420832 -10.9336635
## [41] -8.9248527 -13.7687929 -10.8808798 -14.9080849 -12.0130570
## [46] -13.2631016 -14.1075759 -8.9262800 -13.2798822 -12.0143070
## [51] -19.8939725 -21.6888817 -18.5067252 -12.4549799 -13.9023779
## [56] -15.8533153 -13.9405227 -12.2972153 -16.8346161 -11.2349117
## [61] -14.4207233 -16.7182236 -13.9828581 -18.4659902 -12.4266337
## [66] -17.7015821 -13.4087263 -17.4766913 -12.9127916 -13.2269329
## [71] -12.5587063 -15.0282294 -15.9243220 -16.7978132 -20.0879876
## [76] -15.9938423 -19.8581093 -19.7477287 -13.4007437 -22.2946687
## [81] -15.4248492 -13.9078952 -18.7173113 -13.2457797 -15.5641975
## [86] -13.2318891 -11.7321059 -11.4612747
##
## $dbz
## [1] 11.8458790 11.4754641 10.8560827 9.9854360 8.8621551
## [6] 7.4886835 5.8768395 4.0578211 2.0976111 0.1122841
## [11] -1.7390488 -3.3133638 -4.5663762 -5.5582448 -6.3701775
## [16] -7.0372735 -7.5489517 -7.8837835 -8.0344230 -8.0096686
## [21] -7.8311721 -7.5384172 -7.1927256 -6.8685401 -6.6368757
## [26] -6.5528464 -6.6515867 -6.9494978 -7.4465624 -8.1271327
## [31] -8.9581871 -9.8854455 -10.8302689 -11.6947752 -12.3846543
## [36] -12.8467322 -13.0936848 -13.1902507 -13.2146456 -13.2278985
## [41] -13.2637087 -13.3324027 -13.4305242 -13.5511321 -13.6916541
## [46] -13.8569119 -14.0567183 -14.2998416 -14.5873019 -14.9074910
## [51] -15.2348171 -15.5333097 -15.7659711 -15.9078358 -15.9565105
## [56] -15.9334175 -15.8749606 -15.8195333 -15.7968128 -15.8217940
## [61] -15.8930162 -15.9942147 -16.0995791 -16.1827504 -16.2275006
## [66] -16.2352671 -16.2251647 -16.2268128 -16.2706537 -16.3802168
## [71] -16.5676425 -16.8313863 -17.1546884 -17.5045102 -17.8324754
## [76] -18.0809193 -18.1961576 -18.1453945 -17.9272792 -17.5687318
## [81] -17.1120696 -16.6026318 -16.0823776 -15.5882458 -15.1521143
## [86] -14.8005859 -14.5543805 -14.4276807plot(x = seq(1,len_us), y = newcases_us$positiveIncrease, type = "l", main = 'United States Covid-19 Daily Positive Cases', xlab = 'Time (Days From Start of Pandemic)', ylab = 'Positive Daily Cases')
##Goal Two: Univariate Analysis
###Model Building for Cases in Florida
Stationarity vs non-stationarity - concerns about the data
plotts.sample.wge(newcases_fl$positiveIncrease, arlimits = TRUE)
## $autplt
## [1] 1.00000000 0.90351878 0.85112908 0.81870330 0.78426726 0.72368910
## [7] 0.70892248 0.68040351 0.66167307 0.59052464 0.55736019 0.54019962
## [13] 0.50346164 0.45903473 0.44018864 0.39175866 0.35750776 0.31032227
## [19] 0.29436806 0.25707296 0.21075421 0.16689463 0.15038715 0.11609674
## [25] 0.09827282 0.08510643
##
## $freq
## [1] 0.007407407 0.014814815 0.022222222 0.029629630 0.037037037
## [6] 0.044444444 0.051851852 0.059259259 0.066666667 0.074074074
## [11] 0.081481481 0.088888889 0.096296296 0.103703704 0.111111111
## [16] 0.118518519 0.125925926 0.133333333 0.140740741 0.148148148
## [21] 0.155555556 0.162962963 0.170370370 0.177777778 0.185185185
## [26] 0.192592593 0.200000000 0.207407407 0.214814815 0.222222222
## [31] 0.229629630 0.237037037 0.244444444 0.251851852 0.259259259
## [36] 0.266666667 0.274074074 0.281481481 0.288888889 0.296296296
## [41] 0.303703704 0.311111111 0.318518519 0.325925926 0.333333333
## [46] 0.340740741 0.348148148 0.355555556 0.362962963 0.370370370
## [51] 0.377777778 0.385185185 0.392592593 0.400000000 0.407407407
## [56] 0.414814815 0.422222222 0.429629630 0.437037037 0.444444444
## [61] 0.451851852 0.459259259 0.466666667 0.474074074 0.481481481
## [66] 0.488888889 0.496296296
##
## $db
## [1] 13.9443470 12.1738932 9.4960087 5.0020029 1.2010108
## [6] -1.0130101 -0.1506086 0.2039391 -0.3463438 -0.9726125
## [11] -3.3208756 -5.1446096 -4.6985377 -4.0076175 -3.0331597
## [16] -4.4568722 -6.5691185 -9.5735201 -7.8085456 -4.7203532
## [21] -5.1362477 -6.5033145 -6.5841302 -11.7008911 -15.7768594
## [26] -20.3454088 -14.7914168 -13.3311995 -8.6134078 -11.8506844
## [31] -9.4294479 -10.6796170 -10.8717772 -6.8765678 -5.9591468
## [36] -4.1063148 -6.4365181 -5.2994271 -8.7718362 -14.6806581
## [41] -10.7718204 -10.0417172 -11.3897236 -10.1069047 -12.8544662
## [46] -10.2709895 -10.2140102 -8.6959976 -7.5565160 -7.8110494
## [51] -10.8377249 -12.1847651 -18.4078010 -14.4929054 -22.1053741
## [56] -21.8671829 -18.6457477 -17.1312143 -15.2219426 -8.8827236
## [61] -18.0198672 -11.1478112 -18.8429668 -14.8396876 -10.8187051
## [66] -6.0417196 -6.9676002
##
## $dbz
## [1] 10.8952916 10.4208269 9.6290499 8.5210982 7.1053841
## [6] 5.4084900 3.4951614 1.4961088 -0.3815372 -1.9206026
## [11] -3.0320075 -3.7937247 -4.3311799 -4.7225622 -5.0088115
## [16] -5.2333565 -5.4517814 -5.7160206 -6.0608486 -6.5041876
## [21] -7.0550944 -7.7189474 -8.4927699 -9.3490346 -10.2124871
## [26] -10.9466178 -11.3848655 -11.4222535 -11.0879684 -10.5084533
## [31] -9.8212142 -9.1329494 -8.5201783 -8.0386115 -7.7282189
## [36] -7.6148214 -7.7102817 -8.0113599 -8.4961143 -9.1170132
## [41] -9.7927830 -10.4087161 -10.8436066 -11.0271353 -10.9840189
## [46] -10.8151075 -10.6393942 -10.5538767 -10.6249324 -10.8939499
## [51] -11.3823560 -12.0894051 -12.9791522 -13.9554605 -14.8397014
## [56] -15.4053463 -15.5118058 -15.2066656 -14.6399148 -13.9276436
## [61] -13.1243083 -12.2660250 -11.4056316 -10.6130044 -9.9570803
## [66] -9.4915816 -9.2507536Model IDing of stationary models
Model Building
####Florida Cases - ARMA(2,1) Model
To model Florida cases we start with a base model. We can see that the most favored model by BIC is an AR(1). So we begin by building that model. Our AR(1) has a phi of .975, quite close to the unit circle, which we expect to model strongly wandering behavior.
This model has an AIC estimate of 14.01.
In order to estimate an average ASE, we are running this model over segments of our data with 54 iterations. In each case training with at least seventy data points and predicting on twelve. Overall this produces an average ASE of 6,353,070.
Below our ASE estimates we forecast short term and long term forecasts and both follow AR(1) behavior of data dampening towards our mean.
aic5.wge(newcases_fl$positiveIncrease)## ---------WORKING... PLEASE WAIT...
##
##
## Error in aic calculation at 1 1
## Error in aic calculation at 1 2
## Error in aic calculation at 2 0
## Error in aic calculation at 2 2
## Error in aic calculation at 3 0
## Error in aic calculation at 3 1
## Error in aic calculation at 3 2
## Error in aic calculation at 4 0
## Error in aic calculation at 4 1
## Error in aic calculation at 4 2
## Error in aic calculation at 5 0
## Error in aic calculation at 5 1
## Error in aic calculation at 5 2
## Five Smallest Values of aic## p q aic
## 8 2 1 13.97724
## 4 1 0 14.01245
## 3 0 2 14.80561
## 2 0 1 15.43958
## 1 0 0 16.27428aic5.wge(newcases_fl$positiveIncrease, type = 'bic')## ---------WORKING... PLEASE WAIT...
##
##
## Error in aic calculation at 1 1
## Error in aic calculation at 1 2
## Error in aic calculation at 2 0
## Error in aic calculation at 2 2
## Error in aic calculation at 3 0
## Error in aic calculation at 3 1
## Error in aic calculation at 3 2
## Error in aic calculation at 4 0
## Error in aic calculation at 4 1
## Error in aic calculation at 4 2
## Error in aic calculation at 5 0
## Error in aic calculation at 5 1
## Error in aic calculation at 5 2
## Five Smallest Values of bic## p q bic
## 4 1 0 14.05549
## 8 2 1 14.06332
## 3 0 2 14.87017
## 2 0 1 15.48262
## 1 0 0 16.29580fl_arma_21 = est.arma.wge(newcases_fl$positiveIncrease, p = 2, q = 1)##
## Coefficients of Original polynomial:
## -0.0070 0.9634
##
## Factor Roots Abs Recip System Freq
## 1+0.9850B -1.0152 0.9850 0.5000
## 1-0.9780B 1.0225 0.9780 0.0000
##
## fl_arma_21## $phi
## [1] -0.006958978 0.963353207
##
## $theta
## [1] -0.9221269
##
## $res
## [1] -96.980030 -31.446941 -72.949698 -36.471073 -65.330728
## [6] -45.549146 -61.877488 -36.893978 -64.876273 -32.814665
## [11] -48.465346 -35.457746 -70.039076 -29.808364 29.719937
## [16] -95.862182 -6.379709 -30.423757 -26.217494 131.433577
## [21] -145.518665 -24.636130 363.431958 -282.548857 511.065905
## [26] -478.293879 607.985877 -254.397575 -426.424493 517.214198
## [31] 409.019162 41.458684 -606.888552 167.792148 172.553403
## [36] -673.223311 61.876216 332.217433 -312.910066 -69.985375
## [41] 286.821460 -421.258504 235.449816 -663.078588 632.056247
## [46] 101.816955 -638.099434 47.615189 -6.640751 84.271350
## [51] -455.440684 879.447680 -742.453664 -16.552237 -112.568365
## [56] 50.300458 -384.098389 69.551687 541.043082 -337.470012
## [61] -170.652071 170.581555 -304.057150 -43.831702 246.202377
## [66] -494.646068 362.982017 -193.534920 -306.001397 550.666659
## [71] -496.086986 260.276677 -797.827951 1456.312778 -656.760412
## [76] -179.301368 -177.240334 -230.768626 834.865558 -595.157039
## [81] -31.663234 -75.656780 200.691549 -484.823999 -119.172558
## [86] 171.186415 591.567341 -339.205324 -211.339228 -128.000624
## [91] -74.240039 645.199359 136.823849 -182.026824 -22.061377
## [96] -149.991125 -212.932603 60.319997 290.408665 281.915507
## [101] 231.100730 643.361538 -493.594169 -361.225326 1084.212942
## [106] -165.986553 595.215148 679.105163 257.905827 -499.577450
## [111] -591.625121 383.962847 2259.036748 -308.342849 3850.117969
## [116] 1174.313729 -1202.465388 -2901.551155 485.845142 982.375685
## [121] 3277.089110 111.972557 1580.239818 -560.737458 -4217.026051
## [126] 1487.359042 2462.791359 -446.245537 2181.741834 -281.838325
## [131] 4515.972026 -1516.162835 -4161.355810 1818.911795 3399.515909
##
## $avar
## [1] 1107900
##
## $aic
## [1] 13.97724
##
## $aicc
## [1] 14.9955
##
## $bic
## [1] 14.06332
##
## $se.phi
## [1] 0.0009161559 0.0008599616
##
## $se.theta
## [1] 0.001601026fl_arma_21$aic## [1] 13.97724trainingSize = 70
horizon = 12
ASEHolder = numeric()
for( i in 1:(135-(trainingSize + horizon) + 1))
{
forecasts = fore.aruma.wge(newcases_fl$positiveIncrease[i:(i+(trainingSize-1))],phi = fl_arma_21$phi, theta = fl_arma_21$theta, n.ahead = horizon, plot = FALSE)
ASE = mean((newcases_fl$positiveIncrease[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
ASEHolder[i] = ASE
}
ASEHolder## [1] 167627.68 206397.80 137941.85 555602.91 521952.39
## [6] 102153.75 87478.32 108624.36 183158.65 270619.99
## [11] 71089.92 115747.03 149317.00 133958.81 312340.25
## [16] 400954.76 247387.67 93284.93 158738.58 355742.35
## [21] 730532.31 931261.68 316047.95 439025.90 684184.64
## [26] 1021243.81 1771168.82 2899014.83 3188439.68 2658403.85
## [31] 2249266.64 2963186.04 2424907.70 7652197.48 14157083.00
## [36] 12795601.42 14103207.41 12339977.68 10590941.21 13337047.71
## [41] 19029966.10 29257992.40 30450023.78 15124656.10 17509801.44
## [46] 4594811.59 4654248.02 6537591.44 20214870.57 25595997.81
## [51] 23181259.41 8136658.95 8319256.74 8111410.44#Distribution of ASEs on Two Week Periods
hist(ASEHolder, xlab = "ASE of model at a given Training Set", main = "ASE Distribution for Model ARIMA(2,1,1) for Florida Data")
#Mean ASE
WindowedASE = mean(ASEHolder)
WindowedASE## [1] 6154656##6154656
short_fl_ar1 = fore.aruma.wge(newcases_fl$positiveIncrease,phi = fl_arma_21$phi,theta = fl_arma_21$theta, n.ahead = 7, plot = TRUE)
long_fl_ar1 = fore.aruma.wge(newcases_fl$positiveIncrease,phi = fl_arma_21$phi,theta = fl_arma_21$theta, n.ahead = 90, plot = TRUE)
final_pred = fore.aruma.wge(newcases_fl$positiveIncrease[1:123],phi = fl_arma_21$phi,theta = fl_arma_21$theta, n.ahead = 12, plot = TRUE)
final_pred_df = data.frame(t = seq(124:135), final_pred$f)
plot(newcases_fl$positiveIncrease, type = "l", ylab = "Count of New Cases", xlab = "Time", main = "Florida ARMA(2,1) Model Final 12 Predictions")
lines(ts(final_pred$f, start = 124, end = 135), col = "blue")
MLP/RNN
####MLP Model for Florida Cases
trainingSize = 70
horizon = 12
ASEHolder = numeric()
for( i in 1:(135-(trainingSize + horizon) + 1))
{
mlp.fit = mlp(ts(newcases_fl$positiveIncrease[1:trainingSize+i]), hd = 5, comb = "median")
forecasts = forecast(mlp.fit,h = horizon)
ASE = mean((newcases_fl$positiveIncrease[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] -forecasts$mean)^2)
print(c(i,ASE, "from",trainingSize+i,"to",(trainingSize+ i + (horizon) - 1)))
ASEHolder[i] = ASE
}## [1] "1" "136862.403545455" "from"
## [4] "71" "to" "82"
## [1] "2" "127182.137937909" "from"
## [4] "72" "to" "83"
## [1] "3" "202122.416464837" "from"
## [4] "73" "to" "84"
## [1] "4" "806350.157810554" "from"
## [4] "74" "to" "85"
## [1] "5" "62630.1942021227" "from"
## [4] "75" "to" "86"
## [1] "6" "101744.603764582" "from"
## [4] "76" "to" "87"
## [1] "7" "83378.7696203324" "from"
## [4] "77" "to" "88"
## [1] "8" "83538.9749645744" "from"
## [4] "78" "to" "89"
## [1] "9" "118757.139444703" "from"
## [4] "79" "to" "90"
## [1] "10" "72302.6627071115" "from"
## [4] "80" "to" "91"
## [1] "11" "89839.0388780293" "from"
## [4] "81" "to" "92"
## [1] "12" "96016.3114879308" "from"
## [4] "82" "to" "93"
## [1] "13" "108269.856512192" "from"
## [4] "83" "to" "94"
## [1] "14" "178647.51123949" "from" "84"
## [5] "to" "95"
## [1] "15" "232262.071306686" "from"
## [4] "85" "to" "96"
## [1] "16" "261416.306276133" "from"
## [4] "86" "to" "97"
## [1] "17" "176721.541113281" "from"
## [4] "87" "to" "98"
## [1] "18" "104238.879368051" "from"
## [4] "88" "to" "99"
## [1] "19" "132417.745342906" "from"
## [4] "89" "to" "100"
## [1] "20" "283749.911865593" "from"
## [4] "90" "to" "101"
## [1] "21" "1046628.20717614" "from"
## [4] "91" "to" "102"
## [1] "22" "622389.799876898" "from"
## [4] "92" "to" "103"
## [1] "23" "649469.459830606" "from"
## [4] "93" "to" "104"
## [1] "24" "897502.930280496" "from"
## [4] "94" "to" "105"
## [1] "25" "913673.1442122" "from" "95"
## [5] "to" "106"
## [1] "26" "653069.629102962" "from"
## [4] "96" "to" "107"
## [1] "27" "1980609.47472723" "from"
## [4] "97" "to" "108"
## [1] "28" "2685860.63350751" "from"
## [4] "98" "to" "109"
## [1] "29" "2818067.18379393" "from"
## [4] "99" "to" "110"
## [1] "30" "2824676.96010947" "from"
## [4] "100" "to" "111"
## [1] "31" "2170221.26508948" "from"
## [4] "101" "to" "112"
## [1] "32" "1708444.43668801" "from"
## [4] "102" "to" "113"
## [1] "33" "2724425.51132058" "from"
## [4] "103" "to" "114"
## [1] "34" "6547313.05669394" "from"
## [4] "104" "to" "115"
## [1] "35" "9374037.56095986" "from"
## [4] "105" "to" "116"
## [1] "36" "11497559.8000287" "from"
## [4] "106" "to" "117"
## [1] "37" "11427604.733635" "from" "107"
## [5] "to" "118"
## [1] "38" "7322660.45434971" "from"
## [4] "108" "to" "119"
## [1] "39" "6702484.51188899" "from"
## [4] "109" "to" "120"
## [1] "40" "12031802.3407658" "from"
## [4] "110" "to" "121"
## [1] "41" "17482180.8455421" "from"
## [4] "111" "to" "122"
## [1] "42" "21651704.0440156" "from"
## [4] "112" "to" "123"
## [1] "43" "27845200.687839" "from" "113"
## [5] "to" "124"
## [1] "44" "15146173.9514687" "from"
## [4] "114" "to" "125"
## [1] "45" "1646352645.98057" "from"
## [4] "115" "to" "126"
## [1] "46" "1004816730.24929" "from"
## [4] "116" "to" "127"
## [1] "47" "4706669.00649279" "from"
## [4] "117" "to" "128"
## [1] "48" "51602224.9682025" "from"
## [4] "118" "to" "129"
## [1] "49" "19259729.4836194" "from"
## [4] "119" "to" "130"
## [1] "50" "20740353.3966244" "from"
## [4] "120" "to" "131"
## [1] "51" "16007412.48432" "from" "121"
## [5] "to" "132"
## [1] "52" "15385724.6831372" "from"
## [4] "122" "to" "133"
## [1] "53" "7546847.48620946" "from"
## [4] "123" "to" "134"
## [1] "54" "6720518.71318298" "from"
## [4] "124" "to" "135"ASEHolder## [1] 136862.40 127182.14 202122.42 806350.16 62630.19
## [6] 101744.60 83378.77 83538.97 118757.14 72302.66
## [11] 89839.04 96016.31 108269.86 178647.51 232262.07
## [16] 261416.31 176721.54 104238.88 132417.75 283749.91
## [21] 1046628.21 622389.80 649469.46 897502.93 913673.14
## [26] 653069.63 1980609.47 2685860.63 2818067.18 2824676.96
## [31] 2170221.27 1708444.44 2724425.51 6547313.06 9374037.56
## [36] 11497559.80 11427604.73 7322660.45 6702484.51 12031802.34
## [41] 17482180.85 21651704.04 27845200.69 15146173.95 1646352645.98
## [46] 1004816730.25 4706669.01 51602224.97 19259729.48 20740353.40
## [51] 16007412.48 15385724.68 7546847.49 6720518.71#Distribution of ASEs on Two Week Periods
hist(ASEHolder, xlab = "ASE of model at a given Training Set", main = "ASE Distribution for MLP Model Florida Data")
#Mean ASE
WindowedASE = mean(ASEHolder)
WindowedASE## [1] 54913353plot(forecasts)
#Final Forecasts with data known
mlp.fit_fl_final = mlp(ts(newcases_fl$positiveIncrease[1:123]), hd = 5, comb = "median")
forecasts_fl_mlp = forecast(mlp.fit,h = 12)
all_f = c(rep(1,4), forecasts$fitted, forecasts$mean)
plot(newcases_fl$positiveIncrease, type = "l", ylab = "New Daily Cases", xlab = "Day", main = "")
lines(all_f, col = "blue")
mlp.fit_fl_full = mlp(ts(newcases_fl$positiveIncrease), hd = 5)
forecasts_fl_short = forecast(mlp.fit,h = 7)
forecasts_fl_long = forecast(mlp.fit,h = 90)####Florida Cases Ensemble Model
#ASE fits for ensemble
mlp.fit_fl_final = mlp(ts(newcases_fl$positiveIncrease[1:123]), hd = 5, comb = "median")
forecasts_fl_mlp = forecast(mlp.fit,h = 12)
forecasts_fl_arma = fore.aruma.wge(newcases_fl$positiveIncrease[i:(i+(trainingSize-1))],phi = fl_arma_21$phi, theta = fl_arma_21$theta, n.ahead = 12, plot = FALSE)
ensemble_fl_fore = (forecasts_fl_mlp$mean + forecasts_fl_arma$f) / 2
ensemble_ASE = mean((newcases_fl$positiveIncrease[124:135] - ensemble_fl_fore)^2)
ensemble_ASE## [1] 7287892#8.4 Mill
plot(newcases_fl$positiveIncrease, type = "l", ylab = "Count of New Cases", xlab = "Time", main = "Florida Ensemble Model Final 12 Predictions")
lines(ensemble_fl_fore, col = "blue")
Comparing and Assessing Models
###Model Building for Cases United States
Stationarity vs non-stationarity - concerns about the data
Model IDing of stationary modls
####US Cases ARMA(1,2)
aic5.wge(newcases_us$positiveIncrease)## ---------WORKING... PLEASE WAIT...
##
##
## Error in aic calculation at 1 1
## Error in aic calculation at 2 0
## Error in aic calculation at 2 1
## Error in aic calculation at 2 2
## Error in aic calculation at 3 0
## Error in aic calculation at 3 1
## Error in aic calculation at 3 2
## Error in aic calculation at 4 0
## Error in aic calculation at 4 1
## Error in aic calculation at 4 2
## Error in aic calculation at 5 0
## Error in aic calculation at 5 1
## Error in aic calculation at 5 2
## Five Smallest Values of aic## p q aic
## 6 1 2 15.95786
## 4 1 0 15.96008
## 3 0 2 17.57501
## 2 0 1 18.44652
## 1 0 0 19.55796aic5.wge(newcases_us$positiveIncrease, type = 'bic')## ---------WORKING... PLEASE WAIT...
##
##
## Error in aic calculation at 1 1
## Error in aic calculation at 2 0
## Error in aic calculation at 2 1
## Error in aic calculation at 2 2
## Error in aic calculation at 3 0
## Error in aic calculation at 3 1
## Error in aic calculation at 3 2
## Error in aic calculation at 4 0
## Error in aic calculation at 4 1
## Error in aic calculation at 4 2
## Error in aic calculation at 5 0
## Error in aic calculation at 5 1
## Error in aic calculation at 5 2
## Five Smallest Values of bic## p q bic
## 4 1 0 15.99597
## 6 1 2 16.02964
## 3 0 2 17.62884
## 2 0 1 18.48240
## 1 0 0 19.57590us_arma_12 = est.arma.wge(newcases_us$positiveIncrease, p = 1, q = 2)##
## Coefficients of Original polynomial:
## 0.9919
##
## Factor Roots Abs Recip System Freq
## 1-0.9919B 1.0081 0.9919 0.0000
##
## us_arma_12## $phi
## [1] 0.9919401
##
## $theta
## [1] -0.18789321 -0.06864411
##
## $res
## [1] -210.30561 -105.19842 -127.43409 -130.47115 -128.37416
## [6] -128.55969 -128.66878 -127.63555 -129.81414 -128.48378
## [11] -127.58420 -127.83648 -129.84277 -128.45655 -129.57929
## [16] -125.47155 -131.14212 -128.37475 -129.50547 -125.49104
## [21] -131.14352 -128.37315 -129.50567 -125.49111 -125.14350
## [26] -129.45215 -120.66644 -135.95697 -126.73543 -125.41043
## [31] -131.26822 -126.27465 -125.79469 -119.20347 -125.39426
## [36] -131.64320 -130.06023 -114.94482 -146.78879 -60.97207
## [41] -146.38665 -111.28506 -42.87213 -178.43476 -135.93365
## [46] -43.71046 -95.31929 -38.57411 -70.88751 -56.11484
## [51] 13.02075 188.26586 -299.64272 225.11586 144.71745
## [56] 1852.36521 -932.49698 1411.97729 1258.72068 185.34902
## [61] 2146.62614 1679.93024 -1342.41623 2301.96958 4384.74267
## [66] 386.78843 298.19385 -154.24311 2247.29209 2400.18728
## [71] 511.15148 2004.28971 3605.29519 587.84261 -7953.56392
## [76] 4885.36167 1499.63077 -734.02229 4004.44029 -654.24496
## [81] -3826.72528 -1868.39333 -2034.47714 1049.45065 4654.41605
## [86] -239.01058 854.00880 -4049.06688 270.91090 -1309.20594
## [91] 735.28271 2581.31849 2411.90518 1850.19746 1327.59857
## [96] -8769.73135 -3785.31636 4553.73687 1472.72772 1856.53884
## [101] 3064.16341 -4469.38568 -2714.97901 -2381.03492 579.64339
## [106] 2680.83232 1930.81212 -371.88073 -2796.80904 -2716.64248
## [111] -2553.89407 5078.68403 -2227.91628 5591.16475 -2839.00643
## [116] 178.15218 -4331.14236 1623.50212 -139.18146 387.14321
## [121] 5141.00412 -2898.48024 -2737.02150 -759.39272 -976.02781
## [126] -1943.61036 3160.58072 2750.16971 207.83099 40.11459
## [131] -2114.92169 -816.60334 -131.48326 410.17550 461.03180
## [136] 2404.54211 -702.86684 -4291.39907 -966.35329 614.62840
## [141] 3555.64162 524.76130 1118.26545 1645.44863 -4311.32250
## [146] -2155.25520 5613.45629 -441.53835 3382.40793 2974.00485
## [151] 245.61226 -4864.06610 688.33391 6264.67302 4543.16860
## [156] -745.51228 5330.49330 -1391.99789 -1839.11166 -4767.16587
## [161] 9030.29307 7484.31739 188.46366 1351.13580 -4456.45609
## [166] -8486.20002 6910.13848 3840.80225 9440.58925 -5009.21144
## [171] 8414.72935 -4499.70477 -1414.96658 -1608.42153 5122.92381
## [176] 1720.00550 5811.27282
##
## $avar
## [1] 8142984
##
## $aic
## [1] 15.95786
##
## $aicc
## [1] 16.97115
##
## $bic
## [1] 16.02964
##
## $se.phi
## [1] 0.00008963278
##
## $se.theta
## [1] 0.008017396 0.009972818us_arma_12$aic## [1] 15.95786trainingSize = 70
horizon = 12
ASEHolder = numeric()
for( i in 1:(177-(trainingSize + horizon) + 1))
{
forecasts = fore.aruma.wge(newcases_us$positiveIncrease[i:(i+(trainingSize-1))],phi = us_arma_12$phi, theta = us_arma_12$theta, n.ahead = horizon, plot = FALSE)
ASE = mean((newcases_us$positiveIncrease[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
ASEHolder[i] = ASE
}
ASEHolder## [1] 43725973 35579734 15816563 12379608 16077405 63231777 9750780
## [8] 7422003 7827994 34240666 28525306 6458571 16896630 48641664
## [15] 34877228 15300586 16212745 21221178 19457996 22136822 33803340
## [22] 26312425 16649648 31616011 58839891 87361504 13603088 51265548
## [29] 10098585 16291054 32135503 88968495 21834753 7680905 13124772
## [36] 10367473 12061241 30181559 28890359 8390464 11498325 35163223
## [43] 5970811 9211896 41114819 15001253 16565415 11439564 7133678
## [50] 7371227 7024688 47162686 15221181 3704723 5158051 10417502
## [57] 28158001 5428329 12670962 14728236 15172694 4829085 6396597
## [64] 6776219 6560676 6840061 15007220 12262542 35490861 65797401
## [71] 63948556 24692884 30161102 40875386 45953055 147521940 259676867
## [78] 135521118 159048021 109843644 102674520 144975958 355419374 383704325
## [85] 176877267 88357915 117944780 72816056 113223793 210599000 446366735
## [92] 160337015 55348712 62897740 66463336 160673404#Distribution of ASEs on Two Week Periods
hist(ASEHolder, xlab = "ASE of model at a given Training Set", main = "ASE Distribution for Model ARIMA(2,1,1) for Florida Data")
#Mean ASE
WindowedASE = mean(ASEHolder)
WindowedASE## [1] 55171440#55171440
short_us_ar1 = fore.aruma.wge(newcases_us$positiveIncrease,phi = us_arma_12$phi,theta = us_arma_12$theta, n.ahead = 7, plot = TRUE)
long_us_ar1 = fore.aruma.wge(newcases_us$positiveIncrease,phi = us_arma_12$phi,theta = us_arma_12$theta, n.ahead = 90, plot = TRUE)
final_pred = fore.aruma.wge(newcases_us$positiveIncrease[1:165],phi = us_arma_12$phi,theta = us_arma_12$theta, n.ahead = 12, plot = TRUE)
final_pred_df = data.frame(t = seq(166:177), final_pred$f)
plot(newcases_us$positiveIncrease, type = "l", ylab = "Count of New Cases", xlab = "Time", main = "United States ARMA(1,2) Model Final 12 Predictions")
lines(ts(final_pred$f, start = 166, end = 177), col = "blue")
MLP/RNN
mlp.fit = mlp(ts(newcases_us$positiveIncrease), hd.auto.type = "cv")
mlp.fit## MLP fit with 5 hidden nodes and 20 repetitions.
## Univariate lags: (1,3,4)
## Forecast combined using the median operator.
## MSE: 7065250.643.plot(mlp.fit)
#best option is 20 reps and 5 hds, lags at 1,3, and 4####US Cases MLP
trainingSize = 70
horizon = 12
ASEHolder = numeric()
for( i in 1:(177-(trainingSize + horizon) + 1))
{
mlp.fit = mlp(ts(newcases_us$positiveIncrease[1:trainingSize+i]), hd = 5, reps = 20, lags = c(1,3,4), comb = "median")
forecasts = forecast(mlp.fit,h = horizon)
ASE = mean((newcases_us$positiveIncrease[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] -forecasts$mean)^2)
print(c(i,ASE, "from",trainingSize+i,"to",(trainingSize+ i + (horizon) - 1)))
ASEHolder[i] = ASE
}## [1] "1" "89022533.3070869" "from"
## [4] "71" "to" "82"
## [1] "2" "238946404.28733" "from" "72"
## [5] "to" "83"
## [1] "3" "225286836.558741" "from"
## [4] "73" "to" "84"
## [1] "4" "444016454.95866" "from" "74"
## [5] "to" "85"
## [1] "5" "12200210.9583134" "from"
## [4] "75" "to" "86"
## [1] "6" "12417643.5210084" "from"
## [4] "76" "to" "87"
## [1] "7" "11901125.4391884" "from"
## [4] "77" "to" "88"
## [1] "8" "13030000.3790496" "from"
## [4] "78" "to" "89"
## [1] "9" "163968857.242807" "from"
## [4] "79" "to" "90"
## [1] "10" "120895415.757817" "from"
## [4] "80" "to" "91"
## [1] "11" "39831184.6527887" "from"
## [4] "81" "to" "92"
## [1] "12" "8724324.93672219" "from"
## [4] "82" "to" "93"
## [1] "13" "6637585.18338376" "from"
## [4] "83" "to" "94"
## [1] "14" "9782143.36366572" "from"
## [4] "84" "to" "95"
## [1] "15" "9101342.86342669" "from"
## [4] "85" "to" "96"
## [1] "16" "14867388.0342452" "from"
## [4] "86" "to" "97"
## [1] "17" "17830683.0078988" "from"
## [4] "87" "to" "98"
## [1] "18" "15152727.8489557" "from"
## [4] "88" "to" "99"
## [1] "19" "14475764.1443311" "from"
## [4] "89" "to" "100"
## [1] "20" "14682291.8666076" "from"
## [4] "90" "to" "101"
## [1] "21" "14905521.213616" "from" "91"
## [5] "to" "102"
## [1] "22" "15172229.9477525" "from"
## [4] "92" "to" "103"
## [1] "23" "18953864.1536718" "from"
## [4] "93" "to" "104"
## [1] "24" "25052460.6990796" "from"
## [4] "94" "to" "105"
## [1] "25" "33583290.6272953" "from"
## [4] "95" "to" "106"
## [1] "26" "19856743.6806628" "from"
## [4] "96" "to" "107"
## [1] "27" "12958812.8451124" "from"
## [4] "97" "to" "108"
## [1] "28" "14567321.3031952" "from"
## [4] "98" "to" "109"
## [1] "29" "19955629.7481435" "from"
## [4] "99" "to" "110"
## [1] "30" "31565612.8679613" "from"
## [4] "100" "to" "111"
## [1] "31" "38736757.3103165" "from"
## [4] "101" "to" "112"
## [1] "32" "40614333.5817584" "from"
## [4] "102" "to" "113"
## [1] "33" "34535381.3976825" "from"
## [4] "103" "to" "114"
## [1] "34" "27220659.8306586" "from"
## [4] "104" "to" "115"
## [1] "35" "26334986.8631056" "from"
## [4] "105" "to" "116"
## [1] "36" "34916453.2921645" "from"
## [4] "106" "to" "117"
## [1] "37" "39731732.0974038" "from"
## [4] "107" "to" "118"
## [1] "38" "44900213.2450526" "from"
## [4] "108" "to" "119"
## [1] "39" "42869270.7875341" "from"
## [4] "109" "to" "120"
## [1] "40" "34759923.5768865" "from"
## [4] "110" "to" "121"
## [1] "41" "24640002.5715304" "from"
## [4] "111" "to" "122"
## [1] "42" "30914973.9208757" "from"
## [4] "112" "to" "123"
## [1] "43" "29781114.3698327" "from"
## [4] "113" "to" "124"
## [1] "44" "39400480.5707553" "from"
## [4] "114" "to" "125"
## [1] "45" "43866365.6010222" "from"
## [4] "115" "to" "126"
## [1] "46" "47819504.7876529" "from"
## [4] "116" "to" "127"
## [1] "47" "39861556.5397128" "from"
## [4] "117" "to" "128"
## [1] "48" "37624071.1406232" "from"
## [4] "118" "to" "129"
## [1] "49" "34711568.1202967" "from"
## [4] "119" "to" "130"
## [1] "50" "33505574.599059" "from" "120"
## [5] "to" "131"
## [1] "51" "39829195.8522544" "from"
## [4] "121" "to" "132"
## [1] "52" "41334309.8131301" "from"
## [4] "122" "to" "133"
## [1] "53" "34188235.8501218" "from"
## [4] "123" "to" "134"
## [1] "54" "33474612.2259702" "from"
## [4] "124" "to" "135"
## [1] "55" "22998047.1143707" "from"
## [4] "125" "to" "136"
## [1] "56" "11421740.3774417" "from"
## [4] "126" "to" "137"
## [1] "57" "4098716.16617039" "from"
## [4] "127" "to" "138"
## [1] "58" "15355826.8736451" "from"
## [4] "128" "to" "139"
## [1] "59" "30418661.1507649" "from"
## [4] "129" "to" "140"
## [1] "60" "8812974.05374332" "from"
## [4] "130" "to" "141"
## [1] "61" "6324831.79042253" "from"
## [4] "131" "to" "142"
## [1] "62" "5499983.90616096" "from"
## [4] "132" "to" "143"
## [1] "63" "9064549.01048913" "from"
## [4] "133" "to" "144"
## [1] "64" "3582861.95646854" "from"
## [4] "134" "to" "145"
## [1] "65" "8197331.23287743" "from"
## [4] "135" "to" "146"
## [1] "66" "7465081.55692671" "from"
## [4] "136" "to" "147"
## [1] "67" "7019443.59494492" "from"
## [4] "137" "to" "148"
## [1] "68" "8884151.82972753" "from"
## [4] "138" "to" "149"
## [1] "69" "12676051.142669" "from" "139"
## [5] "to" "150"
## [1] "70" "19853545.2902838" "from"
## [4] "140" "to" "151"
## [1] "71" "21536187.3591235" "from"
## [4] "141" "to" "152"
## [1] "72" "22539786.7924728" "from"
## [4] "142" "to" "153"
## [1] "73" "35957757.345363" "from" "143"
## [5] "to" "154"
## [1] "74" "57810071.5811095" "from"
## [4] "144" "to" "155"
## [1] "75" "83760402.5658058" "from"
## [4] "145" "to" "156"
## [1] "76" "124830117.328055" "from"
## [4] "146" "to" "157"
## [1] "77" "168629733.314794" "from"
## [4] "147" "to" "158"
## [1] "78" "218284200.349663" "from"
## [4] "148" "to" "159"
## [1] "79" "227599793.446626" "from"
## [4] "149" "to" "160"
## [1] "80" "107497989.774778" "from"
## [4] "150" "to" "161"
## [1] "81" "207736742.918578" "from"
## [4] "151" "to" "162"
## [1] "82" "353543947.125516" "from"
## [4] "152" "to" "163"
## [1] "83" "335302463.536474" "from"
## [4] "153" "to" "164"
## [1] "84" "476752798.485755" "from"
## [4] "154" "to" "165"
## [1] "85" "462912417.140798" "from"
## [4] "155" "to" "166"
## [1] "86" "459127801.324589" "from"
## [4] "156" "to" "167"
## [1] "87" "95960053.1530749" "from"
## [4] "157" "to" "168"
## [1] "88" "224759577.860908" "from"
## [4] "158" "to" "169"
## [1] "89" "559958564.656459" "from"
## [4] "159" "to" "170"
## [1] "90" "712516391.431325" "from"
## [4] "160" "to" "171"
## [1] "91" "209174201.182435" "from"
## [4] "161" "to" "172"
## [1] "92" "78634576.4689175" "from"
## [4] "162" "to" "173"
## [1] "93" "86609356.753465" "from" "163"
## [5] "to" "174"
## [1] "94" "63040598.2477316" "from"
## [4] "164" "to" "175"
## [1] "95" "73669051.8690151" "from"
## [4] "165" "to" "176"
## [1] "96" "203055731.319052" "from"
## [4] "166" "to" "177"ASEHolder## [1] 89022533 238946404 225286837 444016455 12200211 12417644 11901125
## [8] 13030000 163968857 120895416 39831185 8724325 6637585 9782143
## [15] 9101343 14867388 17830683 15152728 14475764 14682292 14905521
## [22] 15172230 18953864 25052461 33583291 19856744 12958813 14567321
## [29] 19955630 31565613 38736757 40614334 34535381 27220660 26334987
## [36] 34916453 39731732 44900213 42869271 34759924 24640003 30914974
## [43] 29781114 39400481 43866366 47819505 39861557 37624071 34711568
## [50] 33505575 39829196 41334310 34188236 33474612 22998047 11421740
## [57] 4098716 15355827 30418661 8812974 6324832 5499984 9064549
## [64] 3582862 8197331 7465082 7019444 8884152 12676051 19853545
## [71] 21536187 22539787 35957757 57810072 83760403 124830117 168629733
## [78] 218284200 227599793 107497990 207736743 353543947 335302464 476752798
## [85] 462912417 459127801 95960053 224759578 559958565 712516391 209174201
## [92] 78634576 86609357 63040598 73669052 203055731#Distribution of ASEs on Two Week Periods
hist(ASEHolder, xlab = "ASE of model at a given Training Set", main = "ASE Distribution for MLP Model Florida Data")
#Mean ASE
WindowedASE = mean(ASEHolder)
WindowedASE## [1] 87685290#228 mill
plot(forecasts)
#Actual Forecasting on last segment of data
mlp.fit = mlp(ts(newcases_us$positiveIncrease[1:165]), hd = 5, comb = "median")
forecasts_us_mlp = forecast(mlp.fit,h = 12)
ASE = mean((newcases_us$positiveIncrease[166:177] -forecasts$mean)^2)
ASE## [1] 203055731#53,843,551
plot(newcases_us$positiveIncrease, type = "l")
lines(forecasts$fitted, col = "blue")
all_f = c(forecasts_us_mlp$fitted, forecasts$mean)
plot(newcases_us$positiveIncrease, type = "l")
lines(all_f, col = "blue")
Ensemble
#ASE fits for ensemble
mlp.fit_us_final = mlp(ts(newcases_us$positiveIncrease[1:165]), hd = 5, comb = "median")
forecasts_us_mlp = forecast(mlp.fit_us_final,h = 12)
forecasts_arma_us = fore.aruma.wge(newcases_us$positiveIncrease[1:165],phi = us_arma_12$phi,theta = us_arma_12$theta, n.ahead = 12, plot = TRUE)
ensemble_fore = (forecasts_us_mlp$mean + forecasts_arma_us$f) / 2
ensemble_ASE = mean((newcases_us$positiveIncrease[166:177] -ensemble_fore)^2)
ensemble_ASE## [1] 94399125#114 mill
plot(newcases_us$positiveIncrease, type = "l", ylab = "Count of New Cases", xlab = "Time")
lines(ensemble_fore, col = "blue")
Model Building
##Goal Three: Multivariate Analysis
newcases_fl_multi = initial_data_fl %>% dplyr::select(positiveIncrease, totalTestResultsIncrease, hospitalizedIncrease, deathIncrease)
#Forecast beyond data for Florida
#Forecast future variables
fit.mlp.1 = mlp(ts(newcases_fl_multi$totalTestResultsIncrease),reps = 20, comb = "median")
plot(fit.mlp.1)
fore.mlp.1 = forecast(fit.mlp.1, h = 100)
plot(fore.mlp.1)
fit.mlp.2 = mlp(ts(newcases_fl_multi$hospitalizedIncrease),reps = 20, comb = "median")
plot(fit.mlp.2)
fore.mlp.2 = forecast(fit.mlp.2, h = 100)
plot(fore.mlp.2)
fit.mlp.3 = mlp(ts(newcases_fl_multi$deathIncrease),reps = 20, comb = "median")
plot(fit.mlp.3)
fore.mlp.3 = forecast(fit.mlp.3, h = 100)
plot(fore.mlp.3)
#package them up in data frame.
newvar_fore_fl = data.frame(totalTestResultsIncrease = ts(c(newcases_fl_multi$totalTestResultsIncrease,fore.mlp.1$mean)), hospitalizedIncrease = ts(c(newcases_fl_multi$hospitalizedIncrease,fore.mlp.2$mean)), deathIncrease = ts(c(newcases_fl_multi$deathIncrease,fore.mlp.3$mean)))
#Data has 100 instances beyond current data
dim(newvar_fore_fl)## [1] 235 3###Multivariate Model Building for Florida Cases
####Florida MLR Model
fit = lm(positiveIncrease~totalTestResultsIncrease + hospitalizedIncrease, data = newcases_fl_multi)
summary(fit)##
## Call:
## lm(formula = positiveIncrease ~ totalTestResultsIncrease + hospitalizedIncrease,
## data = newcases_fl_multi)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6794.8 -1047.8 1.2 1357.6 3712.4
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -1358.579010 250.514502 -5.423
## totalTestResultsIncrease 0.137289 0.009719 14.126
## hospitalizedIncrease 5.588560 1.626045 3.437
## Pr(>|t|)
## (Intercept) 0.00000027 ***
## totalTestResultsIncrease < 0.0000000000000002 ***
## hospitalizedIncrease 0.000787 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1687 on 132 degrees of freedom
## Multiple R-squared: 0.7585, Adjusted R-squared: 0.7548
## F-statistic: 207.3 on 2 and 132 DF, p-value: < 0.00000000000000022est_tests = mean(tail(newcases_fl_multi$totalTestResultsIncrease))
est_hospital= mean(tail(newcases_fl_multi$hospitalizedIncrease))
newdata = data.frame(totalTestResultsIncrease = rep(est_tests,12), hospitalizedIncrease = rep(est_hospital,12))
preds = predict(fit, newdata = newdata)
aic5.wge(fit$residuals)#picks 1,1## ---------WORKING... PLEASE WAIT...
##
##
## Five Smallest Values of aic## p q aic
## 5 1 1 13.87980
## 8 2 1 13.89429
## 6 1 2 13.89497
## 11 3 1 13.90310
## 9 2 2 13.90630est1 = est.arma.wge(fit$residuals, p = 1, q = 1)##
## Coefficients of Original polynomial:
## 0.9626
##
## Factor Roots Abs Recip System Freq
## 1-0.9626B 1.0388 0.9626 0.0000
##
## forecasts = fore.arma.wge(fit$residuals,phi = est1$phi,theta = est1$theta, lastn = FALSE,n.ahead = 12)
FinalPredictions_fl_MLR = preds + forecasts$f
plot(newcases_fl$positiveIncrease, type = "l", main ="Florida Multivariate MLR with Correlated Residuals Model", xlab = "Time", ylab = "Count of Cases")
lines(ts(FinalPredictions_fl_MLR, start = 124), col = "blue")
ASE = mean((newcases_fl_multi$positiveIncrease[124:135] - FinalPredictions_fl_MLR)^2)
ASE## [1] 7223923#7223923####Florida Multivariate MLP Cases Model
newcases_fl_multi = initial_data_fl %>% dplyr::select(positiveIncrease, totalTestResultsIncrease, hospitalizedIncrease, deathIncrease)
newcases_fl_var = cbind(ts(newcases_fl_multi$totalTestResultsIncrease),ts(newcases_fl_multi$hospitalizedIncrease),ts(newcases_fl_multi$deathIncrease))
trainingSize = 70
horizon = 12
ASEHolder = numeric()
#Out of bounds if it goes for 54 runs, this ASE will be slightly less wide than the others. But the windowed portion means its average ASE is still good
for( i in 1:(135-(trainingSize + horizon) ))
{
mlp.fit = mlp(ts(newcases_fl_multi$positiveIncrease[1:trainingSize+i]), hd = 5, comb = "median", xreg = newcases_fl_var[1:trainingSize+i,])
forecasts = forecast(mlp.fit,h = horizon, xreg = newcases_fl_var[1:(trainingSize + i + 12),])
ASE = mean((newcases_fl_multi$positiveIncrease[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] -forecasts$mean)^2)
print(c(i,ASE, "from",trainingSize+i,"to",(trainingSize+ i + (horizon) - 1)))
ASEHolder[i] = ASE
}## [1] "1" "1261483.79644511" "from"
## [4] "71" "to" "82"
## [1] "2" "273696.554721272" "from"
## [4] "72" "to" "83"
## [1] "3" "579875.799528204" "from"
## [4] "73" "to" "84"
## [1] "4" "468976.054053524" "from"
## [4] "74" "to" "85"
## [1] "5" "124513.288720101" "from"
## [4] "75" "to" "86"
## [1] "6" "302874.212697658" "from"
## [4] "76" "to" "87"
## [1] "7" "428357.033618399" "from"
## [4] "77" "to" "88"
## [1] "8" "1324153.7769922" "from" "78"
## [5] "to" "89"
## [1] "9" "604107.352111642" "from"
## [4] "79" "to" "90"
## [1] "10" "71017.5449695825" "from"
## [4] "80" "to" "91"
## [1] "11" "114361.143692688" "from"
## [4] "81" "to" "92"
## [1] "12" "140652.239185186" "from"
## [4] "82" "to" "93"
## [1] "13" "268968.434484137" "from"
## [4] "83" "to" "94"
## [1] "14" "237699.720458527" "from"
## [4] "84" "to" "95"
## [1] "15" "286696.723393135" "from"
## [4] "85" "to" "96"
## [1] "16" "448164.572180672" "from"
## [4] "86" "to" "97"
## [1] "17" "161148.484525536" "from"
## [4] "87" "to" "98"
## [1] "18" "231282.062815193" "from"
## [4] "88" "to" "99"
## [1] "19" "185108.116614573" "from"
## [4] "89" "to" "100"
## [1] "20" "368991.631123234" "from"
## [4] "90" "to" "101"
## [1] "21" "489468.02926218" "from" "91"
## [5] "to" "102"
## [1] "22" "720856.803557143" "from"
## [4] "92" "to" "103"
## [1] "23" "616007.958322201" "from"
## [4] "93" "to" "104"
## [1] "24" "1080277.06785036" "from"
## [4] "94" "to" "105"
## [1] "25" "1122319.77989348" "from"
## [4] "95" "to" "106"
## [1] "26" "925910.211210712" "from"
## [4] "96" "to" "107"
## [1] "27" "2233987.80008569" "from"
## [4] "97" "to" "108"
## [1] "28" "2950822.34832411" "from"
## [4] "98" "to" "109"
## [1] "29" "2750720.09652519" "from"
## [4] "99" "to" "110"
## [1] "30" "2977914.31209935" "from"
## [4] "100" "to" "111"
## [1] "31" "2032503.77683334" "from"
## [4] "101" "to" "112"
## [1] "32" "1683654.75681333" "from"
## [4] "102" "to" "113"
## [1] "33" "2957762.82869165" "from"
## [4] "103" "to" "114"
## [1] "34" "6180853.87216433" "from"
## [4] "104" "to" "115"
## [1] "35" "10824963.3463907" "from"
## [4] "105" "to" "116"
## [1] "36" "13311046.6426884" "from"
## [4] "106" "to" "117"
## [1] "37" "9584539.26148209" "from"
## [4] "107" "to" "118"
## [1] "38" "8597984.02861702" "from"
## [4] "108" "to" "119"
## [1] "39" "8532558.4120635" "from" "109"
## [5] "to" "120"
## [1] "40" "9040044.51526702" "from"
## [4] "110" "to" "121"
## [1] "41" "15645700.9258271" "from"
## [4] "111" "to" "122"
## [1] "42" "20174799.4879275" "from"
## [4] "112" "to" "123"
## [1] "43" "15873230.937122" "from" "113"
## [5] "to" "124"
## [1] "44" "11311186.6901928" "from"
## [4] "114" "to" "125"
## [1] "45" "102952915.413264" "from"
## [4] "115" "to" "126"
## [1] "46" "369622968.423746" "from"
## [4] "116" "to" "127"
## [1] "47" "68194283.1503566" "from"
## [4] "117" "to" "128"
## [1] "48" "25300313.5423627" "from"
## [4] "118" "to" "129"
## [1] "49" "13116074.6741099" "from"
## [4] "119" "to" "130"
## [1] "50" "17920209.2743727" "from"
## [4] "120" "to" "131"
## [1] "51" "19745785.4020462" "from"
## [4] "121" "to" "132"
## [1] "52" "14768358.828979" "from" "122"
## [5] "to" "133"
## [1] "53" "8317827.75800403" "from"
## [4] "123" "to" "134"ASEHolder## [1] 1261483.80 273696.55 579875.80 468976.05 124513.29
## [6] 302874.21 428357.03 1324153.78 604107.35 71017.54
## [11] 114361.14 140652.24 268968.43 237699.72 286696.72
## [16] 448164.57 161148.48 231282.06 185108.12 368991.63
## [21] 489468.03 720856.80 616007.96 1080277.07 1122319.78
## [26] 925910.21 2233987.80 2950822.35 2750720.10 2977914.31
## [31] 2032503.78 1683654.76 2957762.83 6180853.87 10824963.35
## [36] 13311046.64 9584539.26 8597984.03 8532558.41 9040044.52
## [41] 15645700.93 20174799.49 15873230.94 11311186.69 102952915.41
## [46] 369622968.42 68194283.15 25300313.54 13116074.67 17920209.27
## [51] 19745785.40 14768358.83 8317827.76#Distribution of ASEs on Two Week Periods
hist(ASEHolder, xlab = "ASE of model at a given Training Set", main = "ASE Distribution for MLP Model Florida Data")
#Mean ASE
WindowedASE = mean(ASEHolder)
WindowedASE## [1] 15083773#18757436 - 18 mill
plot(forecasts)
#Final Forecasts with data known
mlp.fit = mlp(ts(newcases_fl_multi$positiveIncrease[1:123]), hd = 5, comb = "median", xreg =newcases_fl_var[1:123,] )
forecasts = forecast(mlp.fit,h = 12, xreg = newcases_fl_var[1:135,])
fl_multi_mlp_fore = forecasts$mean
all_f = c(forecasts$fitted, forecasts$mean)
plot(newcases_fl_multi$positiveIncrease, type = "l", main = "Florida Multivariate MLP Model with Fits and Final 12 Predictions", xlab = "Time", ylab = "Count of Cases")
lines(all_f, col = "blue")
#Forecast beyond data####Florida Multivariate Ensemble Model
ensemble_fore = (fl_multi_mlp_fore + FinalPredictions_fl_MLR)/2
plot(newcases_fl_multi$positiveIncrease, type = "l", main = "Florida Multivariate Ensemble Model with Final 12 Predictions", xlab = "Time", ylab = "Count of Cases")
lines(ensemble_fore, col = "blue")
ASE_fl_multi = mean((newcases_fl_multi$positiveIncrease[124:135] -ensemble_fore)^2)
ASE_fl_multi## [1] 7845838#ASE of 8,427,522Compare Multivariate Models
###Multivariate US Models
#Forecast variablesnewcases_us_multi = initial_data_us %>% dplyr::select(positiveIncrease, totalTestResultsIncrease, hospitalizedIncrease, deathIncrease)
#Forecast Future
#Forecast future variables
fit.mlp.1 = mlp(ts(newcases_us_multi$totalTestResultsIncrease),reps = 20, comb = "median")
plot(fit.mlp.1)
fore.mlp.1 = forecast(fit.mlp.1, h = 100)
plot(fore.mlp.1)
fit.mlp.2 = mlp(ts(newcases_us_multi$hospitalizedIncrease),reps = 20, comb = "median")
plot(fit.mlp.2)
fore.mlp.2 = forecast(fit.mlp.2, h = 100)
plot(fore.mlp.2)
fit.mlp.3 = mlp(ts(newcases_us_multi$deathIncrease),reps = 20, comb = "median")
plot(fit.mlp.3)
fore.mlp.3 = forecast(fit.mlp.3, h = 100)
plot(fore.mlp.3)
#package them up in data frame.
newvar_fore_us = data.frame(totalTestResultsIncrease = ts(c(newcases_us_multi$totalTestResultsIncrease,fore.mlp.1$mean)), hospitalizedIncrease = ts(c(newcases_us_multi$hospitalizedIncrease,fore.mlp.2$mean)), deathIncrease = ts(c(newcases_us_multi$deathIncrease,fore.mlp.3$mean)))
dim(newvar_fore_us)## [1] 277 3####US MLR with Correlated Errors Model
fit = lm(positiveIncrease~totalTestResultsIncrease + hospitalizedIncrease, data = newcases_us_multi[1:165,])
summary(fit)##
## Call:
## lm(formula = positiveIncrease ~ totalTestResultsIncrease + hospitalizedIncrease,
## data = newcases_us_multi[1:165, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -42222 -3101 -2878 5316 14897
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 3101.227341 869.662497 3.566
## totalTestResultsIncrease 0.049977 0.002809 17.793
## hospitalizedIncrease 2.336046 0.290613 8.038
## Pr(>|t|)
## (Intercept) 0.000477 ***
## totalTestResultsIncrease < 0.0000000000000002 ***
## hospitalizedIncrease 0.00000000000018 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7349 on 162 degrees of freedom
## Multiple R-squared: 0.743, Adjusted R-squared: 0.7399
## F-statistic: 234.2 on 2 and 162 DF, p-value: < 0.00000000000000022est_tests = mean(tail(newcases_us_multi$totalTestResultsIncrease))
est_hospital= mean(tail(newcases_us_multi$hospitalizedIncrease))
newdata = data.frame(totalTestResultsIncrease = rep(est_tests,12), hospitalizedIncrease = rep(est_hospital,12))
preds = predict(fit, newdata = newdata)
aic5.wge(fit$residuals)#picks 3,2 with full data## ---------WORKING... PLEASE WAIT...
##
##
## Five Smallest Values of aic## p q aic
## 5 1 1 16.83224
## 6 1 2 16.83708
## 11 3 1 16.83910
## 8 2 1 16.83932
## 13 4 0 16.83996est1 = est.arma.wge(fit$residuals, p = 3, q = 2)##
## Coefficients of Original polynomial:
## 1.9444 -0.9016 -0.0473
##
## Factor Roots Abs Recip System Freq
## 1-1.9919B+0.9961B^2 0.9998+-0.0652i 0.9981 0.0104
## 1+0.0474B -21.0791 0.0474 0.5000
##
## forecasts = fore.arma.wge(fit$residuals,phi = est1$phi,theta = est1$theta, lastn = FALSE,n.ahead = 12)
FinalPredictions_us_MLR = preds + forecasts$f
plot(newcases_us$positiveIncrease, type = "l", main ="US Multivariate MLR with Correlated Residuals Model", xlab = "Time", ylab = "Count of Cases")
lines(ts(FinalPredictions_us_MLR, start = 166), col = "blue")
ASE = mean((newcases_us_multi$positiveIncrease[166:177] - FinalPredictions_us_MLR)^2)
ASE## [1] 108181241#108181241####US MLP/RNN Model
trainingSize = 70
horizon = 12
ASEHolder = numeric()
for( i in 1:(177-(trainingSize + horizon) + 1))
{
mlp.fit = mlp(ts(newcases_us_multi$positiveIncrease[1:trainingSize+i]), hd = 5, comb = "median", xreg = newvar_fore_us[1:trainingSize+i,])
forecasts = forecast(mlp.fit,h = horizon, xreg = newvar_fore_us[1:(trainingSize + i + 13),])
ASE = mean((newcases_us_multi$positiveIncrease[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] -forecasts$mean)^2)
print(c(i,ASE, "from",trainingSize+i,"to",(trainingSize+ i + (horizon) - 1)))
ASEHolder[i] = ASE
}## [1] "1" "43000778.6911506" "from"
## [4] "71" "to" "82"
## [1] "2" "204625617.705568" "from"
## [4] "72" "to" "83"
## [1] "3" "1466778018.54207" "from"
## [4] "73" "to" "84"
## [1] "4" "517634058.291322" "from"
## [4] "74" "to" "85"
## [1] "5" "1073428892.6305" "from" "75"
## [5] "to" "86"
## [1] "6" "28040619.7675119" "from"
## [4] "76" "to" "87"
## [1] "7" "34587590.668675" "from" "77"
## [5] "to" "88"
## [1] "8" "132034541.796117" "from"
## [4] "78" "to" "89"
## [1] "9" "251980174.586674" "from"
## [4] "79" "to" "90"
## [1] "10" "203771833.098283" "from"
## [4] "80" "to" "91"
## [1] "11" "190487087.955692" "from"
## [4] "81" "to" "92"
## [1] "12" "51994952.7860718" "from"
## [4] "82" "to" "93"
## [1] "13" "35029154.8138601" "from"
## [4] "83" "to" "94"
## [1] "14" "39372017.7627948" "from"
## [4] "84" "to" "95"
## [1] "15" "12565899.971205" "from" "85"
## [5] "to" "96"
## [1] "16" "16304836.0156536" "from"
## [4] "86" "to" "97"
## [1] "17" "24178893.0408641" "from"
## [4] "87" "to" "98"
## [1] "18" "38208456.3770153" "from"
## [4] "88" "to" "99"
## [1] "19" "38021182.6200925" "from"
## [4] "89" "to" "100"
## [1] "20" "33176379.7822393" "from"
## [4] "90" "to" "101"
## [1] "21" "25942733.1175506" "from"
## [4] "91" "to" "102"
## [1] "22" "12104900.7956612" "from"
## [4] "92" "to" "103"
## [1] "23" "15848731.6739233" "from"
## [4] "93" "to" "104"
## [1] "24" "27763094.1123379" "from"
## [4] "94" "to" "105"
## [1] "25" "27895210.5007287" "from"
## [4] "95" "to" "106"
## [1] "26" "44660066.9113287" "from"
## [4] "96" "to" "107"
## [1] "27" "31502234.205533" "from" "97"
## [5] "to" "108"
## [1] "28" "24769370.0605083" "from"
## [4] "98" "to" "109"
## [1] "29" "25635193.4105226" "from"
## [4] "99" "to" "110"
## [1] "30" "37799064.8860735" "from"
## [4] "100" "to" "111"
## [1] "31" "52800134.0676121" "from"
## [4] "101" "to" "112"
## [1] "32" "73508414.3873967" "from"
## [4] "102" "to" "113"
## [1] "33" "59016137.9679405" "from"
## [4] "103" "to" "114"
## [1] "34" "46143816.6579463" "from"
## [4] "104" "to" "115"
## [1] "35" "47439505.7506728" "from"
## [4] "105" "to" "116"
## [1] "36" "51121076.1621967" "from"
## [4] "106" "to" "117"
## [1] "37" "55274536.2972929" "from"
## [4] "107" "to" "118"
## [1] "38" "64994918.4602247" "from"
## [4] "108" "to" "119"
## [1] "39" "69286695.5256278" "from"
## [4] "109" "to" "120"
## [1] "40" "61735688.5508376" "from"
## [4] "110" "to" "121"
## [1] "41" "49787250.055183" "from" "111"
## [5] "to" "122"
## [1] "42" "51705832.7415391" "from"
## [4] "112" "to" "123"
## [1] "43" "53477102.4154954" "from"
## [4] "113" "to" "124"
## [1] "44" "67682125.780084" "from" "114"
## [5] "to" "125"
## [1] "45" "86699430.2144218" "from"
## [4] "115" "to" "126"
## [1] "46" "96655529.5202114" "from"
## [4] "116" "to" "127"
## [1] "47" "97275023.3110608" "from"
## [4] "117" "to" "128"
## [1] "48" "89466185.2146659" "from"
## [4] "118" "to" "129"
## [1] "49" "79837280.3056706" "from"
## [4] "119" "to" "130"
## [1] "50" "72107403.6619167" "from"
## [4] "120" "to" "131"
## [1] "51" "80532367.4369871" "from"
## [4] "121" "to" "132"
## [1] "52" "90729385.4085078" "from"
## [4] "122" "to" "133"
## [1] "53" "98110913.7649145" "from"
## [4] "123" "to" "134"
## [1] "54" "93296151.7410589" "from"
## [4] "124" "to" "135"
## [1] "55" "87797725.265264" "from" "125"
## [5] "to" "136"
## [1] "56" "65097815.4663319" "from"
## [4] "126" "to" "137"
## [1] "57" "22126749.9807346" "from"
## [4] "127" "to" "138"
## [1] "58" "32317939.3810602" "from"
## [4] "128" "to" "139"
## [1] "59" "96120071.7458489" "from"
## [4] "129" "to" "140"
## [1] "60" "105977617.616046" "from"
## [4] "130" "to" "141"
## [1] "61" "101709635.217617" "from"
## [4] "131" "to" "142"
## [1] "62" "102662805.245706" "from"
## [4] "132" "to" "143"
## [1] "63" "66370838.8542523" "from"
## [4] "133" "to" "144"
## [1] "64" "63752793.5612299" "from"
## [4] "134" "to" "145"
## [1] "65" "86278992.7044238" "from"
## [4] "135" "to" "146"
## [1] "66" "100729029.510123" "from"
## [4] "136" "to" "147"
## [1] "67" "98132975.5477682" "from"
## [4] "137" "to" "148"
## [1] "68" "100860039.505746" "from"
## [4] "138" "to" "149"
## [1] "69" "66907459.0858618" "from"
## [4] "139" "to" "150"
## [1] "70" "37211750.6032156" "from"
## [4] "140" "to" "151"
## [1] "71" "77589003.9726429" "from"
## [4] "141" "to" "152"
## [1] "72" "89077029.2104401" "from"
## [4] "142" "to" "153"
## [1] "73" "72360140.6966313" "from"
## [4] "143" "to" "154"
## [1] "74" "77596531.1690199" "from"
## [4] "144" "to" "155"
## [1] "75" "68766942.5563394" "from"
## [4] "145" "to" "156"
## [1] "76" "51514189.1925356" "from"
## [4] "146" "to" "157"
## [1] "77" "60277832.3932788" "from"
## [4] "147" "to" "158"
## [1] "78" "68742260.1176055" "from"
## [4] "148" "to" "159"
## [1] "79" "77462848.1580835" "from"
## [4] "149" "to" "160"
## [1] "80" "57513604.8681221" "from"
## [4] "150" "to" "161"
## [1] "81" "82501776.3109904" "from"
## [4] "151" "to" "162"
## [1] "82" "93608354.701688" "from" "152"
## [5] "to" "163"
## [1] "83" "94330366.0664351" "from"
## [4] "153" "to" "164"
## [1] "84" "211012657.290596" "from"
## [4] "154" "to" "165"
## [1] "85" "170474837.321788" "from"
## [4] "155" "to" "166"
## [1] "86" "160649937.487974" "from"
## [4] "156" "to" "167"
## [1] "87" "214282002.226501" "from"
## [4] "157" "to" "168"
## [1] "88" "160042149.880698" "from"
## [4] "158" "to" "169"
## [1] "89" "149864496.511027" "from"
## [4] "159" "to" "170"
## [1] "90" "501891827.356067" "from"
## [4] "160" "to" "171"
## [1] "91" "496211786.014006" "from"
## [4] "161" "to" "172"
## [1] "92" "108009247.279602" "from"
## [4] "162" "to" "173"
## [1] "93" "98167665.7203687" "from"
## [4] "163" "to" "174"
## [1] "94" "96676019.930379" "from" "164"
## [5] "to" "175"
## [1] "95" "121120949.225492" "from"
## [4] "165" "to" "176"
## [1] "96" "233279323.230828" "from"
## [4] "166" "to" "177"ASEHolder## [1] 43000779 204625618 1466778019 517634058 1073428893 28040620
## [7] 34587591 132034542 251980175 203771833 190487088 51994953
## [13] 35029155 39372018 12565900 16304836 24178893 38208456
## [19] 38021183 33176380 25942733 12104901 15848732 27763094
## [25] 27895211 44660067 31502234 24769370 25635193 37799065
## [31] 52800134 73508414 59016138 46143817 47439506 51121076
## [37] 55274536 64994918 69286696 61735689 49787250 51705833
## [43] 53477102 67682126 86699430 96655530 97275023 89466185
## [49] 79837280 72107404 80532367 90729385 98110914 93296152
## [55] 87797725 65097815 22126750 32317939 96120072 105977618
## [61] 101709635 102662805 66370839 63752794 86278993 100729030
## [67] 98132976 100860040 66907459 37211751 77589004 89077029
## [73] 72360141 77596531 68766943 51514189 60277832 68742260
## [79] 77462848 57513605 82501776 93608355 94330366 211012657
## [85] 170474837 160649937 214282002 160042150 149864497 501891827
## [91] 496211786 108009247 98167666 96676020 121120949 233279323#Distribution of ASEs on Two Week Periods
hist(ASEHolder, xlab = "ASE of model at a given Training Set", main = "ASE Distribution for MLP Model Florida Data")
#Mean ASE
WindowedASE = mean(ASEHolder)
WindowedASE## [1] 117967734#97494363
plot(forecasts)
#Final Forecasts with data known
mlp.fit = mlp(ts(newcases_us_multi$positiveIncrease[1:177]), hd = 5, comb = "median", xreg = newvar_fore_us[1:177,])
forecasts_us_mlp = forecast(mlp.fit,h = 12, xreg = newvar_fore_us[1:190,])
all_f = c(rep(1,4),forecasts_us_mlp$fitted, forecasts_us_mlp$mean)
plot(newcases_us_multi$positiveIncrease, type = "l")
lines(all_f, col = "blue")
#final 12 forecasts
mlp.fit = mlp(ts(newcases_us_multi$positiveIncrease[1:165]), hd = 5, comb = "median", xreg = newvar_fore_us[1:165,])
forecasts_us_mlp = forecast(mlp.fit,h = 12, xreg = newvar_fore_us[1:177,])
all_f = c(rep(1,4),forecasts_us_mlp$fitted, forecasts_us_mlp$mean)
plot(newcases_us_multi$positiveIncrease, type = "l")
lines(all_f, col = "blue")
ASE_final12 = mean((newcases_us_multi$positiveIncrease[166:177] -forecasts_us_mlp$mean)^2)
ASE_final12## [1] 59623538#45799110####US Ensemble
ensemble_us_fore = (forecasts_us_mlp$mean + FinalPredictions_us_MLR)/2
plot(newcases_us_multi$positiveIncrease, type = "l")
lines(ensemble_us_fore, col = "blue")
#Final 12 ASE
ASE_final12 = mean((newcases_us_multi$positiveIncrease[166:177] -ensemble_us_fore)^2)
ASE_final12## [1] 81111113#70596024##Forecasting with new data
#final_us_data = read.csv()
#final_fl_data = read.csv()